function SAT_SST_Step_01_pair_with_CMIP6(P)

    % *********************************************************************
    % Load GHCNmV4
    % *********************************************************************
    Data = SAT_SST_func_load_station_data(P);

    % *********************************************************************
    % Pair HadSST4
    % *********************************************************************
    Data.HadSST4_anm  = SAT_SST_func_pair_HadSST4(Data.lon,Data.lat,P,0);

    % *********************************************************************
    % Pair CMIP6 observations with land data
    % *********************************************************************
    P.do_CMIP6_sens = 0;
    [Data.CMIP6_AT,Data.CMIP6_SST] = SAT_SST_func_pair_CMIP6(Data.lon,Data.lat,P);

    % *********************************************************************
    % Save data
    % *********************************************************************
    Data.T            = permute(Data.T,          [3 1 2]);
    Data.HadSST4_anm  = permute(Data.HadSST4_anm,[3 1 2]);
    Data.CMIP6_AT     = permute(Data.CMIP6_AT,   [3 1 2]);
    Data.CMIP6_SST    = permute(Data.CMIP6_SST,  [3 1 2]);

    % Subset stations that have HadSST4 paired ----------------------------
    %                         and at least 120 months of observational data
    l_use_1 = ~all(isnan(Data.HadSST4_anm(:,:)),2) & (nansum(~isnan(Data.T(:,:)),2)>=120);
    
    % Subset only high quality stations in further analysis
    l_use_2 = SAT_SST_func_get_high_quality_stations(Data.T,Data.HadSST4_anm,P);
    Data    = struct_subset(Data,l_use_1 & l_use_2,1);
    Data.CMIP6_SST_anm = Data.CMIP6_SST - nanmean(Data.CMIP6_SST(:,:,[1982:2014]-1849),3);
    clear('l_use_1','l_use_2')
    
    % save data -----------------------------------------------------------
    raw_SAT_full  = Data.CMIP6_AT;          raw_SAT_full(isnan(Data.T)) = nan;
    raw_SST_anm   = Data.CMIP6_SST_anm;     raw_SST_anm(isnan(Data.HadSST4_anm)) = nan;
    lon           = Data.lon;
    lat           = Data.lat;
    stations      = Data.stations;
    dis           = Data.dis;
    land_fraction = Data.land_fraction;
    file_save = [SAT_SST_func_IO('data',P),'Raw_SAT_SST_',P.date,'.mat'];
    save(file_save,'raw_SAT_full','raw_SST_anm','lon','lat','stations',...
                       'dis','land_fraction','Data','-v7.3');
end
